#include "multi_pack.h"
int Multi_Pack::compare_sp(SingleP &sp_A, SingleP &sp_B) {
  int val = 0;
  int size_A = sp_A.state.size();
  int size_B = sp_B.state.size();
  if (size_A != size_B) {
    val = -1;
    cout << " Wrong @ Multi_Pack::compare_sp sp_A != sp_B" << endl;
    cout << "size_A = " << size_A << "\t size_B = " << size_B << endl;
    sp_A.print();
    cout << "----------- -----------" << endl;
    sp_B.print();

    exit(0);
  }
  for (int i = 0; i < size_A; i++) {
    if (sp_A.state[i].n != sp_B.state[i].n)
      val = -1;
    if (sp_A.state[i].l != sp_B.state[i].l)
      val = -1;
    if (sp_A.state[i].E != sp_B.state[i].E)
      val = -1;
    if (sp_A.state[i].j_2 != sp_B.state[i].j_2)
      val = -1;
    if (sp_A.state[i].tz_2 != sp_B.state[i].tz_2)
      val = -1;
    if (sp_A.state[i].index != sp_B.state[i].index)
      val = -1;
  }
  if (val < 0) {
    cout << " Wrong @ Multi_Pack::compare_sp sp_A != sp_B 2" << endl;

    sp_A.print();
    cout << "----------- -----------" << endl;
    sp_B.print();
    exit(0);
  }
}

void Multi_Pack::compare_tp(TripleP &tp_T, TripleP &tp_x, vector<int> &ch_2x,
                            vector<vector<int>> &index_2x) {
  int ch_size = tp_T.channel.size();
  ch_2x.resize(ch_size);
  index_2x.resize(ch_size);

  for (int ch = 0; ch < ch_size; ch++) {
    int J_2 = tp_T.channel[ch].J_2;
    int T_2 = tp_T.channel[ch].T_2;
    int ch_2x_t = -1;
    for (int ch_x = 0; ch_x < tp_x.channel.size(); ch_x++) {
      if (tp_x.channel[ch_x].J_2 == J_2 && tp_x.channel[ch_x].T_2 == T_2) {
        ch_2x_t = ch_x;
        break;
      }
    }
    ch_2x[ch] = ch_2x_t;
    int ch_x = ch_2x_t;
    // if (ch_x < 0) {
    //   index_2x[ch].resize(0);
    //   continue;
    // }
    int size = tp_T.state[ch].size();
    int size_x = tp_x.state[ch_x].size();
    index_2x[ch].resize(size);
    for (int i = 0; i < size; i++) {
      int x = -1;
      if (ch_x < 0) {
        index_2x[ch][i] = -1;
        continue;
      }
      for (int j = 0; j < size_x; j++) {
        if (tp_T.state[ch][i].a != tp_x.state[ch_x][j].a)
          continue;
        if (tp_T.state[ch][i].b != tp_x.state[ch_x][j].b)
          continue;
        if (tp_T.state[ch][i].c != tp_x.state[ch_x][j].c)
          continue;
        if (tp_T.state[ch][i].Jab_2 != tp_x.state[ch_x][j].Jab_2)
          continue;
        if (tp_T.state[ch][i].Tab_2 != tp_x.state[ch_x][j].Tab_2)
          continue;
        // if (tp_T.state[ch][i].Tab_2z != tp_x.state[ch_x][j].Tab_2z)
        //   continue;
        // if (tp_T.state[ch][i].J_2 != tp_x.state[ch_x][j].J_2)
        //   continue;
        // if (tp_T.state[ch][i].T_2 != tp_x.state[ch_x][j].T_2)
        //   continue;

        x = j;
        break;
      }
      index_2x[ch][i] = x;
    }
  }
}

void Multi_Pack::init_sub(string file, SMatrix_Ch &vec_t,
                          SMatrix_Ch &vec_flag_t, TripleP &triple_t,
                          SingleP &single_t) {
  LEC_CONST LECs;
  // string file_beg = "serial_dat/" + file + "/";
  string file_beg = file + "/";
  string sp_file = file_beg + "SP_serial.dat";
  string tp_file = file_beg + "TP_serial.dat";
  string jaco_file = file_beg + "jaco_serial.dat";
  string anti_file = file_beg + "anti_serial.dat";
  string jacoNL_file = file_beg + "jacoNL_serial.dat";
  string trans_file = file_beg + "trans_serial.dat";
  // string wave_init_file = file_beg + "Wave_init.Input";
  string cont_file = file_beg + "Cont_serial.dat";
  string ope_file = file_beg + "Ope_serial.dat";
  string tpe_file = file_beg + "Tpe_serial.dat";
  string pack_file = file_beg + "Package_serial.dat";
  // string TP_filename = file_beg + "TP" + ".Input";
  string DF_filename = file_beg + "DF.Input";

  Get_Info get_info;
  get_info.init(DF_filename);
  int output_flag = -1;
  get_info.build(output_flag);

  CG_coeff cg_x;
  SixJ_coeff sixJ_x;
  Binom bin_t;
  // cg_x.init(10);

  // sixJ_x.init(12);
  NinJ_coeff ninJ_x;
  // ninJ_x.init(5);
  Hobra hobra_t(bin_t, cg_x, ninJ_x);
  char u1[256]; // unuse

  Serial_SL ser_t;

  // SingleP single_t;
  // cout << "\t \t sp_file : " << sp_file << endl;
  ser_t.load(single_t, sp_file);
  string file_sp_out = file_beg + "sp.dat";
  single_t.print(file_sp_out);
  // single_t.print();
  // TripleP triple_t;
  // cout << "\t \t tp_file : " << tp_file << endl;
  ser_t.load(triple_t, tp_file);
  string file_tp_out = file_beg + "tp.dat";
  string file_tp_ch_out = file_beg + "tp_ch.dat";

  // cout << "file_tp_out : " << file_tp_out << endl;
  // cout << "file_tp_ch_out : " << file_tp_ch_out << endl;
  triple_t.print(file_tp_out);
  triple_t.print_channel(file_tp_ch_out);

  Jacobi jaco_b;
  // ser_t.load(jaco_b, jaco_file);
  Antisym antisym_t(sixJ_x, ninJ_x, hobra_t, jaco_b);
  // ser_t.load(antisym_t, anti_file);
  Jacobi_NL jaco_nl;
  // ser_t.load(jaco_nl, jacoNL_file);
  TransCoeff trans_t(sixJ_x, ninJ_x, hobra_t, jaco_nl, single_t, triple_t);
  // ser_t.load(trans_t, trans_file);
  // cout << "Begin to load Wave" << endl;
  //   double hbar_omg;
  // vector<int> wave_const(13);
  // std::ifstream WaveInit_file;
  // WaveInit_file.open(wave_init_file);
  // for (int i = 0; i < 13; i++) {
  //
  //   int A = -1;
  //   if (i == 1) {
  //     WaveInit_file >> hbar_omg;
  //   } else {
  //     WaveInit_file >> A;
  //   }
  //   // cout << A << "\t** " << endl;
  //   wave_const[i] = A;
  //   WaveInit_file.getline(u1, 256, '\n');
  // }
  // WaveInit_file.close();
  // int lamb_cut = wave_const[0];
  // //double hbar_omg = wave_const[1];
  // int reg_pow = wave_const[2];
  // int k_mesh = wave_const[3];
  // int k_scale = wave_const[4];
  // int p_mesh = wave_const[5];
  // int Emax_F = wave_const[6];
  // int Emax_I = wave_const[7];

  const int lamb_cut = get_info.wave_info.Lambda_Cut;
  double hbar_omg = get_info.wave_info.hbar_omega;
  int reg_pow = get_info.wave_info.regulator_power;
  int k_mesh = get_info.wave_info.k_mesh;
  int k_scale = get_info.wave_info.k_scale;
  int p_mesh = get_info.wave_info.p_mesh;
  int Emax_F = get_info.wave_info.EF_max;
  int Emax_I = get_info.wave_info.EI_max;
  int n_max = get_info.wave_info.n_max;
  int l_max = get_info.wave_info.l_max;
  int lbarS = get_info.wave_info.lbarS;
  int lbarS_G = get_info.wave_info.lbarS_G;
  int lbarS_P = get_info.wave_info.lbarS_P;

  Wave_Fun wave(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale, p_mesh, Emax_F,
                Emax_I);
  Cont_term cont_x(LECs, wave, sixJ_x, jaco_b);

  // ser_t.load(cont_x, cont_file);
  // cont_x.print("data/cont_term.dat");

  Ope_term ope_x(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);

  // ser_t.load(ope_x, ope_file);
  // ope_x.print("data/ope_term.dat");

  Tpe_term tpe_t(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);

  // ser_t.load(tpe_t, tpe_file);
  // cout << "Begin to load Package" << endl;
  //~~~~~~~~~~~~~~~~ Package  ~~~~~~~~~~~~~~~~~~~//
  Package pack(triple_t, jaco_nl, jaco_b, trans_t, antisym_t, cont_x, ope_x,
               tpe_t);
  // cout<<"Package file : "<<pack_file<<endl;
  ser_t.load(pack, pack_file); /**/
  vec_t.clear();
  vec_flag_t.clear();
  vec_t = pack.vec_S;
  vec_flag_t = pack.vec_flag_S;

  string pack_file_dat = file_beg + "package.dat";
  pack.print(pack_file_dat);

  // ~~~~~~~~~~~~~~~~ Some Constant ~~~~~~~~~~~~~~~~~~//

  // int E_cut_t;
  // std::ifstream TPinput_file;
  // TPinput_file.open(TP_filename);
  //
  // for (int i = 0; i < 6; i++) {
  //
  //   int A = -1;
  //   TPinput_file >> A;
  //   // cout << A << "\t** " << endl;
  //   // abc_vec[i] = A;
  //   TPinput_file.getline(u1, 256, '\n');
  // }
  // int TP_allow_file_num;
  // TPinput_file >> TP_allow_file_num;
  // TPinput_file.getline(u1, 256, '\n');
  // TPinput_file >> E_cut_t;
  // TPinput_file.close();
  int E_max_t = get_info.tp_info.E_max;
  cout << " !!! flag_const = " << flag_const << endl;
  if (flag_const == 1) {
    this->hbar_omega = hbar_omg;
    this->regulator_pow = reg_pow;
    this->Cut_Lambda = lamb_cut;
    this->E_max = E_max_t;
    cout << "\t hbar_omega : " << hbar_omega << endl;
    cout << "\t regulator_pow : " << regulator_pow << endl;
    cout << "\t Cut_Lambda : " << Cut_Lambda << endl;
    cout << "\t E_max : " << E_max << endl;
    flag_const = -1;
  } else {
    int flag = 1;
    if (this->hbar_omega != hbar_omg)
      flag = -1;
    if (this->regulator_pow != reg_pow)
      flag = -1;
    if (this->Cut_Lambda != lamb_cut)
      flag = -1;
    if (this->E_max != E_max_t)
      flag = -1;
    if (flag < 0) {
      cout << "Wrong happaned @ Multi_Pack::init_sub" << endl;
      cout << "hw : " << this->hbar_omega << "\t " << hbar_omg << endl;
      cout << "regulator_pow : " << this->regulator_pow << "\t " << reg_pow
           << endl;
      cout << "Cut_Lambda : " << this->Cut_Lambda << "\t " << lamb_cut << endl;
      exit(0);
    }
  }
  // cout << "Hbar_Omega = " << this->hbar_omega << endl;
  cout << "E_max = " << this->E_max << "\t " << E_max_t << endl;
  int flag = 1;
  if (LECs.LEC_C1 != pack.LEC_C1_P)
    flag = -1;
  if (LECs.LEC_C3 != pack.LEC_C3_P)
    flag = -1;
  if (LECs.LEC_C4 != pack.LEC_C4_P)
    flag = -1;
  if (LECs.LEC_CD != pack.LEC_CD_P)
    flag = -1;
  if (LECs.LEC_CE != pack.LEC_CE_P)
    flag = -1;
  if (LECs.LEC_gA != pack.LEC_gA_P)
    flag = -1;
  if (flag < 0) {
    cout << "Wrong @ Multi_Pack::init_sub " << endl;
    cout << "LECs different !!! " << endl;
    cout << "\t pack.LEC\t LECs" << endl;
    cout << " c1 = " << pack.LEC_C1_P << "\t " << LECs.LEC_C1 << endl;
    cout << " c3 = " << pack.LEC_C3_P << "\t " << LECs.LEC_C3 << endl;
    cout << " c4 = " << pack.LEC_C4_P << "\t " << LECs.LEC_C4 << endl;
    cout << " cD = " << pack.LEC_CD_P << "\t " << LECs.LEC_CD << endl;
    cout << " cE = " << pack.LEC_CE_P << "\t " << LECs.LEC_CE << endl;
    cout << " gA = " << pack.LEC_gA_P << "\t " << LECs.LEC_gA << endl;

    exit(0);
  } else {
    cout << " LECs GOOD ! " << endl;
  }
}

void Multi_Pack::init(string file_sub1, string file_sub2) {

  LEC_CONST LECs;
  this->LEC_CD = LECs.LEC_CD;
  this->LEC_CE = LECs.LEC_CE;
  this->LEC_C1 = LECs.LEC_C1;
  this->LEC_C3 = LECs.LEC_C3;
  this->LEC_C4 = LECs.LEC_C4;
  this->LEC_gA = LECs.LEC_gA;

  this->init_sub(file_sub1, this->vec_A_S, this->vec_flag_A_S, TP_A, SP_A);
  cout << "Finished init_sub A " << endl;
  this->init_sub(file_sub2, this->vec_B_S, this->vec_flag_B_S, TP_B, SP_B);
  cout << "Finished init_sub B " << endl;
  // SingleP single_t;

  // this->SP.build_T("data/Input/SP.Input");
  // this->SP.build_T_E("data/Input/SP_E.Input");
  // const std::string path_sp = "data/Input/SP.Input";
  // const std::string path_sp_E = "data/Input/SP_E.Input";

  Get_Info get_info;
  get_info.init("data/Input/DF.Input");
  get_info.build();

  // bf::path path(path_s);
  // if (bf::exists(path_sp)) {
  //   this->SP.build_T(path_sp);
  //   cout << "Build SP using data/Input/SP.Input " << endl;
  // } else if (bf::exists(path_sp_E)) {
  //   this->SP.build_T_E(path_sp_E);
  //   cout << "Build SP using data/Input/SP_E.Input " << endl;
  // }
  if (get_info.sp_info.Type == "Read") {
    this->SP.build_T(get_info.sp_info.sp_vec);
  } else if (get_info.sp_info.Type == "E_max") {
    this->SP.build_T_E(get_info.sp_info.E_max);
  } else {
    cerr << " Wrong @ build SP " << endl;
    cerr << "\t get_info.sp_info.Type : " << get_info.sp_info.Type << endl;
    exit(0);
  }

  // single_t.build_Tz("data/SP_Tz.Input");
  this->SP.build_Tz();

  this->SP.print("serial_dat/Output_MP/Multi_SP.dat");
  this->SP.print_Tz("serial_dat/Output_MP/Multi_SP_Tz.dat");
  // exit(0);

  // single_t.build_Tz(0,1,1,1);
  // single_t.build_TzM(0,1,1,1);

  // vector<int> abc_vec(6);
  // int TP_allow_file_num;
  // char u1[256]; // unuse
  // std::ifstream TPinput_file;
  // TPinput_file.open("data/Input/TP.Input");
  //
  // for (int i = 0; i < 6; i++) {
  //
  //   int A = -1;
  //   TPinput_file >> A;
  //   // cout << A << "\t** " << endl;
  //   abc_vec[i] = A;
  //   TPinput_file.getline(u1, 256, '\n');
  // }
  // TPinput_file >> TP_allow_file_num;
  // TPinput_file.close();

  // int a_min = abc_vec[0];
  // int a_max = abc_vec[1];
  // int b_min = abc_vec[2];
  // int b_max = abc_vec[3];
  // int c_min = abc_vec[4];
  // int c_max = abc_vec[5];

  // vector<vector<vector<int>>> TP_allow_gather;
  // vector<int> E_cut_gather;
  // vector<int> Same_Particles_gather;
  // vector<vector<int>> E_sign_gather;

  //
  // int TP_allow_file_num = 2;

  // std::ifstream TP_Allow_file;
  // for (int i = 0; i < TP_allow_file_num; i++) {
  //   vector<vector<int>> TP_allow_vec;
  //   TP_allow_vec.resize(2);
  //   for (int j = 0; j < 2; j++) {
  //     TP_allow_vec[j].resize(6);
  //   }
  //   vector<int> E_sign(6);
  //   int Same_Particles_Flag;
  //   int E_cut;
  //   int ii = i;
  //   string filename = "data/Input/TP_allow_" + std::to_string(ii);
  //   filename += ".Input";
  //   TP_Allow_file.open(filename);
  //   TP_Allow_file >> E_cut;
  //   TP_Allow_file.getline(u1, 256, '\n');
  //   TP_Allow_file >> Same_Particles_Flag;
  //   TP_Allow_file.getline(u1, 256, '\n');
  //   for (int i = 0; i < 6; i++) {
  //     int min, max;
  //     TP_Allow_file >> min >> max >> E_sign[i];
  //     // cout << "\t TP allow : " << min << "\t" << max << endl;
  //     // cout << A << "\t** " << endl;
  //     TP_allow_vec[0][i] = min;
  //     TP_allow_vec[1][i] = max;
  //     TP_Allow_file.getline(u1, 256, '\n');
  //   }
  //   TP_Allow_file.close();
  //   TP_allow_gather.push_back(TP_allow_vec);
  //   E_cut_gather.push_back(E_cut);
  //   E_sign_gather.push_back(E_sign);
  //   Same_Particles_gather.push_back(Same_Particles_Flag);
  //   this->SP.print(TP_allow_vec);
  // }

  vector<vector<vector<int>>> TP_allow_gather;
  vector<int> E_cut_gather;
  vector<int> Same_Particles_gather;
  vector<vector<int>> E_sign_gather;

  int TP_allow_file_num = get_info.tp_al_info.size();
  for (int i = 0; i < TP_allow_file_num; i++) {
    vector<vector<int>> TP_allow_vec(2);
    for (int j = 0; j < 2; j++) {
      TP_allow_vec[j].resize(6);
      for (int k = 0; k < 6; k++) {
        TP_allow_vec[j][k] = get_info.tp_al_info[i].tp_al_vec[k][j];
      }
    }
    vector<int> E_sign(6);
    for (int k = 0; k < 6; k++) {
      E_sign[k] = get_info.tp_al_info[i].tp_al_vec[k][2];
    }
    int Same_Particles_Flag = get_info.tp_al_info[i].Same_particle_Flag;
    int E_cut = get_info.tp_al_info[i].E_cut;
    TP_allow_gather.push_back(TP_allow_vec);
    E_cut_gather.push_back(E_cut);
    E_sign_gather.push_back(E_sign);
    Same_Particles_gather.push_back(Same_Particles_Flag);
    this->SP.print(TP_allow_vec);
  }

  // ~~~~~~~~~~~~~~~~~~~~~~~ Build Triple needed ~~~~~~~~~~~~~~~~~~~~~~~//
  // TripleP triple_t;
  // cout << "a_max : " << a_max << "\t b_max : " << b_max
  //      << "\t c_max : " << c_max << endl;
  int file_num = E_sign_gather.size();
  /*for (int k = 0; k < file_num; k++) {
    cout << "\t File : " << k << "\t Ecut : " << E_cut_gather[k]
         << "\t Same particles : " << Same_Particles_gather[k] << endl;
    cout << "\t E_sign :" << endl;
    for (int i = 0; i < 6; i++) {
      cout << E_sign_gather[k][i] << "\t";
    }
    cout << endl;
  }*/
  if (get_info.tp_info.Type == "abc") {
    int a_min = get_info.tp_info.abc_vec[0];
    int a_max = get_info.tp_info.abc_vec[1];
    int b_min = get_info.tp_info.abc_vec[2];
    int b_max = get_info.tp_info.abc_vec[3];
    int c_min = get_info.tp_info.abc_vec[4];
    int c_max = get_info.tp_info.abc_vec[5];
    cout << "Start Build T" << endl;
    // this->TP.build_T(this->SP, a_min, a_max, b_min, b_max, c_min, c_max);
    this->TP.init_Constrain(TP_allow_gather);
    this->TP.build_T_Constrain(this->SP, a_min, a_max, b_min, b_max, c_min,
                               c_max);
    cout << "Finished build T" << endl;
  } else if (get_info.tp_info.Type == "E_max") {
    int Emax_t = get_info.tp_info.E_max;
    // this->TP.build_T(this->SP, Emax_t);
    this->TP.init_Constrain(TP_allow_gather);
    this->TP.build_T_Constrain(this->SP, Emax_t);
  } else {
    cerr << " Wrong @ build TP " << endl;
    cerr << "\t get_info.tp_info.Type : " << get_info.tp_info.Type << endl;
    exit(0);
  }
  // this->TP.build_T(this->SP, a_min, a_max, b_min, b_max, c_min, c_max);
  // this->TP.Select_FI_T(this->SP, TP_allow_gather);

  this->TP.Delet_Sort_T();
  Time_Count TC;
  TC.start();
  // cout << "Start Delet_conf_T " << endl;
  this->TP.Delet_conf_T(this->SP, E_sign_gather, E_cut_gather,
                        Same_Particles_gather);
  // cout << "End Delet_conf_T : " << TC.sec() << "s" << endl;

  TC.start();
  cout << "Start build_conf_T " << endl;
  this->TP.build_conf_T(this->SP, E_sign_gather, E_cut_gather,
                        Same_Particles_gather);
  // cout << "End build_conf_T : " << TC.sec() << "s" << endl;

  // this->TP.print_TPconf_info("data/TPconf.dat", "data/TPconf_E.dat",
  //                            "data/TPconf_E_tot.dat");
  // this->TP.sort();
  // this->TP.print("data/Muli_TP_bef.dat");
  cout << "Start pirnt TP " << endl;
  this->TP.print("serial_dat/Output_MP/Multi_TP.dat");

  // exit(0);
  this->TP.print_channel("serial_dat/Output_MP/Multi_TP_channel.dat");

  if (get_info.tp_info.Type == "abc") {
    int a_min = 2 * get_info.tp_info.abc_vec[0];
    int a_max = 2 * get_info.tp_info.abc_vec[1] + 1;
    int b_min = 2 * get_info.tp_info.abc_vec[2];
    int b_max = 2 * get_info.tp_info.abc_vec[3] + 1;
    int c_min = 2 * get_info.tp_info.abc_vec[4];
    int c_max = 2 * get_info.tp_info.abc_vec[5] + 1;
    cout << "Start build TP_Tz " << endl;
    this->TP.build_Tz_Constrain(this->SP, a_min, a_max, b_min, b_max, c_min, c_max);
    cout << "End build TP_Tz " << endl;
  } else if (get_info.tp_info.Type == "E_max") {
    int Emax_t = get_info.tp_info.E_max;
    cout << "Start build TP_Tz " << endl;
    this->TP.build_Tz_Constrain(this->SP, Emax_t);
    cout << "End build TP_Tz " << endl;
  } else {
    cerr << " Wrong @ build TP_Tz " << endl;
    cerr << "\t get_info.tp_info.Type : " << get_info.tp_info.Type << endl;
    exit(0);
  }
  // this->TP.build_Tz(this->SP, 2 * a_min, 2 * a_max + 1, 2 * b_min,
  //                   2 * b_max + 1, 2 * c_min, 2 * c_max + 1);
  //this->TP.Select_FI_Tz(TP_allow_gather);
  this->TP.Delet_Sort_Tz();
  TC.start();
  // cout << "Start Delet_conf_Tz " << endl;
  this->TP.Delet_conf_Tz(this->SP, E_sign_gather, E_cut_gather,
                         Same_Particles_gather);
  // cout << "End Delet_conf_Tz : " << TC.sec() << "s" << endl;
  TC.start();
  // cout << "Start build_conf_Tz_S " << endl;
  this->TP.build_conf_Tz_S(this->SP, E_sign_gather, E_cut_gather,
                           Same_Particles_gather);
  // cout << "End build_conf_Tz_S : " << TC.sec() << "s" << endl;
  // this->TP.print_Tz("data/Muli_TP_Tz_bef.dat");

  // this->TP.sort_Tz();
  // this->TP.print_TPconf_info_Tz("data/TPconf_Tz.dat",
  // "data/TPconf_E_Tz.dat");

  this->TP.print_Tz("serial_dat/Output_MP/Multi_TP_Tz.dat");

  this->TP.print_channel_Tz("serial_dat/Output_MP/Multi_TP_channel_Tz.dat");

  // cout << "check TP.FI_compare(1,100110) : " << TP.FI_compare(1, 100110)
  //     << endl;
  // cout << "check TP.FI_compare(1,100101) : " << TP.FI_compare(1, 100101)
  //     << endl;
  // ~~~~~~~~~~~~~~~~~~~~~~~ Build Triple needed ~~~~~~~~~~~~~~~~~~~~~~~//
  // cout << "Test 0" << endl;
  /*this->init_sub(file_sub1, this->vec_A_S, this->vec_flag_A_S, TP_A, SP_A);
  this->init_sub(file_sub2, this->vec_B_S, this->vec_flag_B_S, TP_B, SP_B);/**/
  // cout << "Test 1" << endl;
  this->compare_sp(SP, SP_A);
  // cout << "Test 2" << endl;
  this->compare_sp(SP, SP_B);

  this->compare_tp(this->TP, this->TP_A, this->ch_2A, this->index_2A);
  this->compare_tp(this->TP, this->TP_B, this->ch_2B, this->index_2B);
  cout << "Finished init Multi_Pack" << endl;
  // this->print_2AB("data/check_multi_pack_ch.dat",
  //                 "data/check_multi_pack_index.dat");
  // cout << "check & print_2AB" << endl;
}

void Multi_Pack::build() {
  LEC_CONST LECs;
  int size_ch = TP.channel.size();
  vec.resize(size_ch);
  vec_flag.resize(size_ch);
  for (int ch = 0; ch < size_ch; ch++) {
    int size = TP.state[ch].size();
    vec[ch].resize(size);
    vec_flag[ch].resize(size);
    // cout << "\t ch : " << ch << "\t size = " << size << endl;
    for (int i = 0; i < size; i++) {
      vec[ch][i].resize(size);
      vec_flag[ch][i].resize(size);
      for (int j = 0; j < size; j++) {
        vec[ch][i][j] = 0.0;
        vec_flag[ch][i][j] = 0;
      }
    }
  }
  for (int ch = 0; ch < size_ch; ch++) {
    int ch_A = ch_2A[ch];
    int ch_B = ch_2B[ch];
    if (ch_A < 0 && ch_B < 0) {
      continue;
    }
    int size = TP.state[ch].size();
    // this->vec_S.cal_beg(size, size);
    // this->vec_flag_S.cal_beg(size, size);
    for (int i = 0; i < size; i++) {
      // this->vec_S.cal_row_off(ch, i);
      // this->vec_flag_S.cal_row_off(ch, i);
      int i_A = -1;
      int i_B = -1;
      if (ch_A >= 0)
        i_A = index_2A[ch][i];
      if (ch_B >= 0)
        i_B = index_2B[ch][i];
      if (i_A < 0 && i_B < 0)
        continue;

      for (int j = 0; j < size; j++) {
        int j_A = -1;
        int j_B = -1;
        if (ch_A >= 0)
          j_A = index_2A[ch][j];
        if (ch_B >= 0)
          j_B = index_2B[ch][j];
        if (j_A < 0 && j_B < 0)
          continue;

        if ((i_A < 0 || j_A < 0) && (i_B < 0 || j_B < 0))
          continue;

        double val = 0.0;
        double val_A = vec_A_S.val(ch_A, i_A, j_A);
        double val_B = vec_B_S.val(ch_B, i_B, j_B);
        int flag_A = vec_flag_A_S.val(ch_A, i_A, j_A);
        int flag_B = vec_flag_B_S.val(ch_B, i_B, j_B);

        if (i_B < 0 || j_B < 0) {
          // vec_A exist
          if (flag_A == 0)
            continue;
          val = val_A;
        } else if (i_A < 0 || j_A < 0) {
          // vec_B exist
          if (flag_B == 0)
            continue;
          val = val_B;
        } else {
          // vec_A & vec_B exist
          if (flag_A == 0 && flag_B == 0)
            continue;

          if (flag_A == 1 && flag_B == 0)
            val = val_A;

          if (flag_A == 0 && flag_B == 1)
            val = val_B;

          if (flag_A == 1 && flag_B == 1) {
            val = val_A;

            if (std::abs(val - val_B) > this->EPS) {
              cout << "Wrong @ Multi_Pack::build 1 " << endl;
              cout << "ch_A : " << ch_A << "\t i_A : " << i_A
                   << "\t j_A : " << j_A << "\t val" << val_A << endl;
              cout << "ch_B : " << ch_B << "\t i_B : " << i_B
                   << "\t j_B : " << j_B << "\t val" << val_B << endl;
              exit(0);
            }
          }
        }
        if (vec_flag[ch][i][j] != 0 &&
            std::abs(vec[ch][i][j] - val) > this->EPS) {
          cout << "Wrong @ Multi_Pack::build 2 " << endl;
          cout << "ch_A : " << ch_A << "\t i_A : " << i_A << "\t j_A : " << j_A
               << "\t val" << vec_A_S.val(ch_A, i_A, j_A) << endl;
          cout << "ch_B : " << ch_B << "\t i_B : " << i_B << "\t j_B : " << j_B
               << "\t val" << vec_B_S.val(ch_B, i_B, j_B) << endl;
          exit(0);
        } /**/
        vec[ch][i][j] = val;
        vec_flag[ch][i][j] = 1;
        // this->vec_S.cal_col_ind(ch, i, j, val);
        // this->vec_flag_S.cal_col_ind(ch, i, j, 1);
      }
    }
  }
  cout << "Finished build Multi_Pack 2" << endl;
  this->check("data/check_multi_pack.dat");
  vec_A_S.clear();
  // vec_flag_A.clear();
  vec_B_S.clear();
  // vec_flag_B.clear();
  // ~~~~~~~~~~~~~~~~~~~~~~ Package build Tz ~~~~~~~~~~~~~~~~~~~~~~~~~~ //
  CG_coeff cg_x;
  SixJ_coeff sixJ_x;
  Binom bin_t;
  cg_x.init(10);

  sixJ_x.init(12);
  NinJ_coeff ninJ_x;
  ninJ_x.init(5, sixJ_x);
  Hobra hobra_t(bin_t, cg_x, ninJ_x);
  char u1[256]; // unuse

  Serial_SL ser_t;

  Jacobi jaco_b;
  // ser_t.load(jaco_b, jaco_file);
  Antisym antisym_t(sixJ_x, ninJ_x, hobra_t, jaco_b);
  // ser_t.load(antisym_t, anti_file);
  Jacobi_NL jaco_nl;
  // ser_t.load(jaco_nl, jacoNL_file);
  TransCoeff trans_t(sixJ_x, ninJ_x, hobra_t, jaco_nl, this->SP, this->TP);
  // ser_t.load(trans_t, trans_file);

  vector<int> wave_const(13);
  std::ifstream WaveInit_file;
  WaveInit_file.open("data/Input/Wave_init.Input");
  for (int i = 0; i < 13; i++) {

    int A = -1;
    WaveInit_file >> A;
    // cout << A << "\t** " << endl;
    wave_const[i] = A;
    WaveInit_file.getline(u1, 256, '\n');
  }
  WaveInit_file.close();
  int lamb_cut = wave_const[0];
  double hbar_omg = wave_const[1];
  int reg_pow = wave_const[2];
  int k_mesh = wave_const[3];
  int k_scale = wave_const[4];
  int p_mesh = wave_const[5];
  int Emax_F = wave_const[6];
  int Emax_I = wave_const[7];
  Wave_Fun wave(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale, p_mesh, Emax_F,
                Emax_I);
  Cont_term cont_x(LECs, wave, sixJ_x, jaco_b);

  // ser_t.load(cont_x, cont_file);
  // cont_x.print("data/cont_term.dat");

  Ope_term ope_x(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);

  // ser_t.load(ope_x, ope_file);
  // ope_x.print("data/ope_term.dat");

  Tpe_term tpe_t(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);

  // ser_t.load(tpe_t, tpe_file);

  //~~~~~~~~~~~~~~~~ Package  ~~~~~~~~~~~~~~~~~~~//
  Package pack(this->TP, jaco_nl, jaco_b, trans_t, antisym_t, cont_x, ope_x,
               tpe_t);
  cout << "Start package.build_S_Tz" << endl;
  vec_S.build_ch(vec);
  vec_flag_S.build_ch(vec_flag);
  vec.clear();
  vec_flag.clear();
  pack.build_S_Tz(cg_x, this->SP, this->vec_S, this->vec_flag_S, this->vec_Tz_S,
                  this->vec_flag_Tz_S);
  // cout << "Start to Rebuild" << endl;

  cout << "Finshed package.build_Tz" << endl;
  // Check Vs nnlosat
  /*
  int a_min = 0;
  int a_max = 9;
  int b_min = 0;
  int b_max = 9;
  int c_min = 0;
  int c_max = 9;
  std::vector<int> abc_range = {a_min, a_max, b_min, b_max, c_min, c_max};
  pack.Readnnlosat_SP(this->SP, "data/sat/Navr_ord_states.dat");
  pack.Compare_satVSzhuo("data/sat/3BME_part_3.dat", "data/sat/satVsYZ.dat",
                         abc_range);/**/
}

void Multi_Pack::build_S() {
  LEC_CONST LECs;
  cout << "Start to Build Multi_Pack" << endl;
  int size_ch = TP.channel.size();
  this->vec_S.cal_init(size_ch, 1);
  this->vec_flag_S.cal_init(size_ch, 0);
  for (int ch = 0; ch < size_ch; ch++) {
    // cout << "ch : " << ch << endl;
    int ch_A = ch_2A[ch];
    int ch_B = ch_2B[ch];
    if (ch_A < 0 && ch_B < 0) {
      continue;
    }
    int size = TP.state[ch].size();
    // cout << "\t \t size = " << size << endl;
    // cout << "\t \t \t index_2A[ch].size() = " << index_2A[ch].size() << endl;
    // cout << "\t \t \t index_2B[ch].size() = " << index_2B[ch].size() << endl;
    this->vec_S.cal_beg(size, size);
    this->vec_flag_S.cal_beg(size, size);
    for (int i = 0; i < size; i++) {
      this->vec_S.cal_row_off(ch, i);
      this->vec_flag_S.cal_row_off(ch, i);
      int i_A = -1;
      int i_B = -1;
      if (ch_A >= 0)
        i_A = index_2A[ch][i];
      if (ch_B >= 0)
        i_B = index_2B[ch][i];
      if (i_A < 0 && i_B < 0)
        continue;

      for (int j = 0; j < size; j++) {
        int j_A = -1;
        int j_B = -1;
        if (ch_A >= 0)
          j_A = index_2A[ch][j];
        if (ch_B >= 0)
          j_B = index_2B[ch][j];
        if (j_A < 0 && j_B < 0)
          continue;

        if ((i_A < 0 || j_A < 0) && (i_B < 0 || j_B < 0))
          continue;

        double val = 0.0;

        if (i_B < 0 || j_B < 0) {
          // vec_A exist
          double val_A = vec_A_S.val(ch_A, i_A, j_A);
          int flag_A = vec_flag_A_S.val(ch_A, i_A, j_A);

          if (flag_A == 0)
            continue;
          val = val_A;
        } else if (i_A < 0 || j_A < 0) {
          // vec_B exist
          double val_B = vec_B_S.val(ch_B, i_B, j_B);
          int flag_B = vec_flag_B_S.val(ch_B, i_B, j_B);

          if (flag_B == 0)
            continue;
          val = val_B;
        } else {
          // vec_A & vec_B exist
          double val_A = vec_A_S.val(ch_A, i_A, j_A);
          int flag_A = vec_flag_A_S.val(ch_A, i_A, j_A);
          double val_B = vec_B_S.val(ch_B, i_B, j_B);
          int flag_B = vec_flag_B_S.val(ch_B, i_B, j_B);
          if (flag_A == 0 && flag_B == 0)
            continue;

          if (flag_A == 1 && flag_B == 0)
            val = val_A;

          if (flag_A == 0 && flag_B == 1)
            val = val_B;

          if (flag_A == 1 && flag_B == 1) {
            val = val_A;

            if (std::abs(val - val_B) > this->EPS) {
              cout << "Wrong @ Multi_Pack::build 1 " << endl;
              cout << "ch_A : " << ch_A << "\t i_A : " << i_A
                   << "\t j_A : " << j_A << "\t val" << val_A << endl;
              cout << "ch_B : " << ch_B << "\t i_B : " << i_B
                   << "\t j_B : " << j_B << "\t val" << val_B << endl;
              exit(0);
            } /**/
          }
        }
        /*if (vec_flag[ch][i][j] != 0 &&
            std::abs(vec[ch][i][j] - val) > this->EPS) {
          cout << "Wrong @ Multi_Pack::build 2 " << endl;
          cout << "ch_A : " << ch_A << "\t i_A : " << i_A << "\t j_A : " << j_A
               << "\t val" << vec_A_S.val(ch_A, i_A, j_A) << endl;
          cout << "ch_B : " << ch_B << "\t i_B : " << i_B << "\t j_B : " << j_B
               << "\t val" << vec_B[ch_B][i_B][j_B] << endl;
          exit(0);
        }/**/
        // vec[ch][i][j] = val;
        // vec_flag[ch][i][j] = 1;
        this->vec_S.cal_col_ind(ch, i, j, val);
        this->vec_flag_S.cal_col_ind(ch, i, j, 1);
      }
    }
    this->vec_S.cal_end(ch);
    this->vec_flag_S.cal_end(ch);
  }
  // cout << "Finished build Multi_Pack 2" << endl;
  // this->check("data/check_multi_pack.dat");
  vec_A_S.clear();
  // vec_flag_A.clear();
  vec_B_S.clear();
  // vec_flag_B.clear();
  // ~~~~~~~~~~~~~~~~~~~~~~ Package build Tz ~~~~~~~~~~~~~~~~~~~~~~~~~~ //
  CG_coeff cg_x;
  SixJ_coeff sixJ_x;
  Binom bin_t;
  cg_x.init(10);

  sixJ_x.init(12);
  NinJ_coeff ninJ_x;
  ninJ_x.init(5, sixJ_x);
  Hobra hobra_t(bin_t, cg_x, ninJ_x);
  char u1[256]; // unuse

  Serial_SL ser_t;

  Jacobi jaco_b;
  // ser_t.load(jaco_b, jaco_file);
  Antisym antisym_t(sixJ_x, ninJ_x, hobra_t, jaco_b);
  // ser_t.load(antisym_t, anti_file);
  Jacobi_NL jaco_nl;
  // ser_t.load(jaco_nl, jacoNL_file);
  TransCoeff trans_t(sixJ_x, ninJ_x, hobra_t, jaco_nl, this->SP, this->TP);
  // ser_t.load(trans_t, trans_file);

  vector<int> wave_const(13);
  std::ifstream WaveInit_file;
  WaveInit_file.open("data/Input/Wave_init.Input");
  for (int i = 0; i < 13; i++) {

    int A = -1;
    WaveInit_file >> A;
    // cout << A << "\t** " << endl;
    wave_const[i] = A;
    WaveInit_file.getline(u1, 256, '\n');
  }
  WaveInit_file.close();
  int lamb_cut = wave_const[0];
  double hbar_omg = wave_const[1];
  int reg_pow = wave_const[2];
  int k_mesh = wave_const[3];
  int k_scale = wave_const[4];
  int p_mesh = wave_const[5];
  int Emax_F = wave_const[6];
  int Emax_I = wave_const[7];
  Wave_Fun wave(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale, p_mesh, Emax_F,
                Emax_I);
  Cont_term cont_x(LECs, wave, sixJ_x, jaco_b);

  // ser_t.load(cont_x, cont_file);
  // cont_x.print("data/cont_term.dat");

  Ope_term ope_x(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);

  // ser_t.load(ope_x, ope_file);
  // ope_x.print("data/ope_term.dat");

  Tpe_term tpe_t(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);

  // ser_t.load(tpe_t, tpe_file);

  //~~~~~~~~~~~~~~~~ Package  ~~~~~~~~~~~~~~~~~~~//
  Package pack(this->TP, jaco_nl, jaco_b, trans_t, antisym_t, cont_x, ope_x,
               tpe_t);
  cout << "Start package.build_S_Tz" << endl;
  pack.build_S_Tz(cg_x, this->SP, this->vec_S, this->vec_flag_S, this->vec_Tz_S,
                  this->vec_flag_Tz_S);
  // cout << "Start to Rebuild" << endl;

  cout << "Finshed package.build_Tz 2" << endl;
  // Check Vs nnlosat
  /*
  int a_min = 0;
  int a_max = 9;
  int b_min = 0;
  int b_max = 9;
  int c_min = 0;
  int c_max = 9;
  std::vector<int> abc_range = {a_min, a_max, b_min, b_max, c_min, c_max};
  pack.Readnnlosat_SP(this->SP, "data/sat/Navr_ord_states.dat");
  pack.Compare_satVSzhuo("data/sat/3BME_part_3.dat", "data/sat/satVsYZ.dat",
                         abc_range);/**/
}

void Multi_Pack::check(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size_ch = vec_S.size_ch;
  int flag = 1;
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_S.col_num_ch[ch];
    int ch_A = ch_2A[ch];
    int ch_B = ch_2B[ch];
    file << "Channel : " << ch << "\t ch_A : " << ch_A << "\t ch_B : " << ch_B
         << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        double val = vec_S.val(ch, i, j);
        if (this->vec_flag_S.val(ch, i, j) == 1) {
          // if (std::abs(val) > this->EPS) {
          file << i << "\t" << j << "\t" << vec_flag_S.val(ch, i, j) << "\t"
               << std::setiosflags(std::ios::fixed) << std::setprecision(3)
               << vec_S.val(ch, i, j) << "\t";

          file << index_2A[ch][i] << "\t" << index_2A[ch][j] << "\t";
          if (index_2A[ch][i] < 0 || index_2A[ch][j] < 0) {
            file << "nan \t";
          } else {
            if (vec_flag_A_S.val(ch_A, index_2A[ch][i], index_2A[ch][j]) == 0) {
              file << "not_cal \t";
            } else {
              double val_A =
                  vec_A_S.val(ch_A, index_2A[ch][i], index_2A[ch][j]);
              file << std::setiosflags(std::ios::fixed) << std::setprecision(3)
                   << val_A << "\t";
              if (std::abs(val_A - val) > this->EPS) {
                flag = -1;
              }
            }
          }

          file << index_2B[ch][i] << "\t" << index_2B[ch][j] << "\t";
          if (index_2B[ch][i] < 0 || index_2B[ch][j] < 0) {
            file << "nan \t";
          } else {
            if (vec_flag_B_S.val(ch_B, index_2B[ch][i], index_2B[ch][j]) == 0) {
              file << "not_cal \t";
            } else {
              double val_B =
                  vec_B_S.val(ch_B, index_2B[ch][i], index_2B[ch][j]);
              file << std::setiosflags(std::ios::fixed) << std::setprecision(3)
                   << val_B << "\t";
              if (std::abs(val_B - val) > this->EPS) {
                flag = -1;
              }
            }
          }

          file << endl;
        }
      }
    }
  }
  file.close();
  if (flag > 0) {
    cout << "Check multi_pack result : Good !" << endl;
  } else {
    cout << "Check multi_pack result : Bad !" << endl;
    exit(0);
  }
}

void Multi_Pack::print_2AB(std::string filename1, std::string filename2) {
  std::ofstream file_ch;
  file_ch.open(filename1);
  int size_ch = ch_2A.size();

  for (int ch = 0; ch < size_ch; ch++) {
    file_ch << "ch : " << ch << "\t ch_A : " << ch_2A[ch]
            << "\t ch_B : " << ch_2B[ch] << endl;
  }
  file_ch.close();

  std::ofstream file_ind;
  file_ind.open(filename2);
  // int size_ch = vec.size();
  for (int ch = 0; ch < size_ch; ch++) {
    file_ind << " --- Channel : " << ch << "\t ch_A : " << ch_2A[ch]
             << "\t ch_B :" << ch_2B[ch] << " --- " << endl;
    int size = index_2A[ch].size();
    for (int i = 0; i < size; i++) {
      file_ind << "\t i : " << i << " \t i_A : " << index_2A[ch][i]
               << " \t i_B : " << index_2B[ch][i] << endl;
    }
  }
  file_ind.close();
}

void Multi_Pack::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size_ch = vec_S.size_ch;
  for (int ch = 0; ch < size_ch; ch++) {
    int size = vec_S.row_num_ch[ch];
    file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (vec_flag_S.val(ch, i, j) == 1)
          file << i << "\t" << j << "\t" << vec_S.val(ch, i, j) << "\t "
               << vec_flag_S.val(ch, i, j) << endl;
      }
    }
  }
  file.close();
}

void Multi_Pack::print_Tz(std::string filename) {
  std::ofstream file;
  file.open(filename);
  // vector<vector<vector<double>>> vec_Tz;
  // vector<vector<vector<int>>> vec_flag_Tz;
  // this->vec_Tz_S.Rebuild(vec_Tz);
  // this->vec_flag_Tz_S.Rebuild(vec_flag_Tz);
  int size_ch = vec_Tz_S.size_ch;

  for (int ch = 0; ch < size_ch; ch++) {
    // cout << "\t ttt " << endl;
    int size = vec_Tz_S.row_num_ch[ch];
    // cout << "Channel : " << ch << "\t size : " << size;
    file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        int a = 0;
        if (vec_flag_Tz_S.val(ch, i, j) == 1)
          file << i << "\t" << j << "\t" << vec_Tz_S.val(ch, i, j) << "\t"
               << vec_flag_Tz_S.val(ch, i, j) << endl;
      }
    }
  }
  file.close();
}

void Multi_Pack::print_General() {
  cout << "  ==========================================================="
       << endl;
  cout << "      __        __           ___  __   __   ___  __   ___    "
       << endl;
  cout << "     |  \\  /\\  |__) |__/    |__  |__) /  \\ |__  /__`  |      "
       << endl;
  cout << "     |__/ /--\\ |  \\ |  \\    |    |  \\ \\__/ |___ .__/  |      "
       << endl;
  cout << "                                                             "
       << endl;
  cout << "                         version 2.0                         "
       << endl;
  cout << "                     by Yuanzhuo 2020-June                   "
       << endl;
  cout << "    ------------------------------------------------------   "
       << endl;
  cout << endl;
  Time_Count TC;
  TC.print_time();
  cout << endl;
}

void Multi_Pack::print_General(std::ofstream &file) {
  file << "  ==========================================================="
       << endl;
  file << "      __        __           ___  __   __   ___  __   ___    "
       << endl;
  file << "     |  \\  /\\  |__) |__/    |__  |__) /  \\ |__  /__`  |     "
       << endl;
  file << "     |__/ /--\\ |  \\ |  \\    |    |  \\ \\__/ |___ .__/  |    "
       << endl;
  file << "                                                             "
       << endl;
  file << "    ------------------------------------------------------   "
       << endl;
  file << endl;
  Time_Count TC;
  TC.print_time(file);
  file << endl;
}

void Multi_Pack::print_TP_ALL_f(std::ofstream &file) {
  Get_Info get_info;
  get_info.init("data/Input/DF.Input");
  get_info.build(-1);

  LEC_CONST LECs;
  file << " ---- Parameters : ----" << endl;
  file << "\t Lambda_Cut = " << get_info.wave_info.Lambda_Cut << "  (MeV)"
       << endl;
  file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  file << "\t cD = " << LECs.LEC_CD << endl;
  file << "\t cE = " << LECs.LEC_CE << endl;
  file << "\t gA = " << LECs.LEC_gA << endl;
  file << "\t Hbar-Omega = " << get_info.wave_info.hbar_omega << endl;
  file << "\t The Power of Regulator : " << get_info.wave_info.regulator_power
       << endl;

  file << " ---- 3BME range ----  " << endl;
  int TP_allow_file_num = get_info.tp_al_info.size();
  for (int i = 0; i < TP_allow_file_num; i++) {
    vector<vector<int>> TP_allow_vec(2);
    for (int j = 0; j < 2; j++) {
      TP_allow_vec[j].resize(6);
      for (int k = 0; k < 6; k++) {
        TP_allow_vec[j][k] = get_info.tp_al_info[i].tp_al_vec[k][j];
      }
    }
    SP.print(TP_allow_vec, file);
  }

  double EPS = 1E-8;
  file << " Only output abs(ME) > " << EPS << endl;
  file << "a\tb\tc\td\te\t f\t J2ab\t J2cd\t J2\t T2ab\t T2cd\t T2\t\t "
          "<abc|V|def>"
       << endl;

  int ch_size = TP.channel.size();
  for (int ch = 0; ch < ch_size; ch++) {
    int size = vec_S.row_num_ch[ch];
    // cout << " \t " << size << "  \t" << TP.state[ch].size() << endl;
    // file << "Channel : " << ch << endl;
    for (int i = 0; i < size; i++) {
      int a = TP.state[ch][i].a;
      int b = TP.state[ch][i].b;
      int c = TP.state[ch][i].c;
      int J2ab = TP.state[ch][i].Jab_2;
      int T2ab = TP.state[ch][i].Tab_2;
      int J2 = TP.state[ch][i].J_2;
      int T2 = TP.state[ch][i].T_2;
      for (int j = 0; j < size; j++) {
        int d = TP.state[ch][j].a;
        int e = TP.state[ch][j].b;
        int f = TP.state[ch][j].c;
        int J2de = TP.state[ch][j].Jab_2;
        int T2de = TP.state[ch][j].Tab_2;
        if (vec_flag_S.val(ch, i, j) == 1) {
          // file << i << "\t" << j << "\t" << vec_S.val(ch, i, j) << "\t "
          //     << vec_flag_S.val(ch, i, j) << endl;
          double val = vec_S.val(ch, i, j);
          if (std::abs(val) < EPS)
            continue;
          file << a << "\t " << b << "\t " << c << "\t ";
          file << d << "\t " << e << "\t " << f << "\t ";
          file << J2ab << "\t " << J2de << "\t " << J2 << "\t ";
          file << T2ab << "\t " << T2de << "\t " << T2 << "\t\t ";
          file << val << endl;
        }
      }
    }
  }
} /**/
